home *** CD-ROM | disk | FTP | other *** search
/ Ham Radio 2000 #1 / Ham Radio 2000.iso / ham2000 / misc / daubwave / daubwave.c next >
Encoding:
C/C++ Source or Header  |  1993-01-15  |  25.5 KB  |  1,001 lines
  1. /* daubwave - This program decomposes a data set by orthogonal wavelets.
  2.           The results are then filtered by setting to zero a number
  3.           of iterations from the low end or from the high end.
  4.           The data is then reconstructed from the modified data.    */
  5.  
  6. #include<stdio.h>
  7. #include<stdlib.h>
  8. /*     stdlib.h is not necessary on all platforms - remove if it causes
  9.     a compiler error and try again                    */
  10. #include<string.h>
  11.  
  12. #define MAXNUM 4096
  13.  
  14. void VectorCopy( float *, float *, long, int );
  15. void Wavelet( float *, long, int, double *, double * );
  16. void VectorPermute( float *, float *, long, int );
  17. void VectorShift( float *, long, int, int );
  18. FILE *FileOpen( char *, char * );
  19. long FileLoad( float *, struct opt );
  20. void FileSave( float *, int, long, struct opt );
  21. void *MemAlloc( int, long );
  22. void TranCoef( double *, double *, int );
  23. long SpanData( int *, long );
  24. char HyperHex( int );
  25. void HelpFile( void );
  26.  
  27. struct opt {   char infile[81];
  28.            char outfile[81];
  29.            char root[81];
  30.            char extension[5];
  31.            int filter;
  32.            int compress;
  33.            char action;
  34.            int level;
  35.            int shift;
  36.            double normal;
  37.            int unix;
  38.            int rotate; };
  39.  
  40. struct opt InitOption( int, char *[] );
  41.  
  42. double coef[11][20] = {      {    1.000000000000, 0.000000000000,
  43.                 0.000000000000, 0.000000000000,
  44.                 0.000000000000, 0.000000000000,
  45.                 0.000000000000, 0.000000000000,
  46.                 0.000000000000, 0.000000000000,
  47.                 0.000000000000, 0.000000000000,
  48.                 0.000000000000, 0.000000000000,
  49.                 0.000000000000, 0.000000000000,
  50.                 0.000000000000, 0.000000000000,
  51.                 0.000000000000, 0.000000000000  },
  52.  
  53.                  {    0.707106781187, 0.707106781187,
  54.                 0.000000000000, 0.000000000000,
  55.                 0.000000000000, 0.000000000000,
  56.                 0.000000000000, 0.000000000000,
  57.                 0.000000000000, 0.000000000000,
  58.                 0.000000000000, 0.000000000000,
  59.                 0.000000000000, 0.000000000000,
  60.                 0.000000000000, 0.000000000000,
  61.                 0.000000000000, 0.000000000000,
  62.                 0.000000000000, 0.000000000000  },
  63.  
  64.                 {    0.482962913145, 0.836516303738,
  65.                 0.224143868042, -0.129409522551,
  66.                 0.000000000000, 0.000000000000,
  67.                 0.000000000000, 0.000000000000,
  68.                 0.000000000000, 0.000000000000,
  69.                 0.000000000000, 0.000000000000,
  70.                 0.000000000000, 0.000000000000,
  71.                 0.000000000000, 0.000000000000,
  72.                 0.000000000000, 0.000000000000,
  73.                 0.000000000000, 0.000000000000  },
  74.  
  75.                 {    0.332670552950, 0.806891509311,
  76.                 0.459877502118, -0.135011020010,
  77.                 -0.085441273882, 0.035226291882,
  78.                 0.000000000000, 0.000000000000,
  79.                 0.000000000000, 0.000000000000,
  80.                 0.000000000000, 0.000000000000,
  81.                 0.000000000000, 0.000000000000,
  82.                 0.000000000000, 0.000000000000,
  83.                 0.000000000000, 0.000000000000,
  84.                 0.000000000000, 0.000000000000  },
  85.  
  86.                 {    0.230377813309, 0.714846570553,
  87.                 0.630880767930, -0.027983769417,
  88.                 -0.187034811719, 0.030841381836,
  89.                 0.032883011667, -0.010597401785,
  90.                 0.000000000000, 0.000000000000,
  91.                 0.000000000000, 0.000000000000,
  92.                 0.000000000000, 0.000000000000,
  93.                 0.000000000000, 0.000000000000,
  94.                 0.000000000000, 0.000000000000,
  95.                 0.000000000000, 0.000000000000  },
  96.  
  97.                 {    0.160102397974, 0.603829269797,
  98.                 0.724308528438, 0.138428145901,
  99.                 -0.242294887066, -0.032244869585,
  100.                 0.077571493840, -0.006241490213,
  101.                 -0.012580751999, 0.003335725285,
  102.                 0.000000000000, 0.000000000000,
  103.                 0.000000000000, 0.000000000000,
  104.                 0.000000000000, 0.000000000000,
  105.                 0.000000000000, 0.000000000000,
  106.                 0.000000000000, 0.000000000000  },
  107.  
  108.                 {    0.111540743350, 0.494623890398,
  109.                 0.751133908021, 0.315250351709,
  110.                 -0.226264693965, -0.129766867567,
  111.                 0.097501605587, 0.027522865530,
  112.                 -0.031582039318, 0.000553842201,
  113.                 0.004777257511, -0.001077301085,
  114.                 0.000000000000, 0.000000000000,
  115.                 0.000000000000, 0.000000000000,
  116.                 0.000000000000, 0.000000000000,
  117.                 0.000000000000, 0.000000000000  },
  118.  
  119.                 {    0.077852054085, 0.396539319482,
  120.                 0.729132090846, 0.469782287405,
  121.                 -0.143906003929, -0.224036184994,
  122.                 0.071309219267, 0.080612609151,
  123.                 -0.038029936935, -0.016574541631,
  124.                 0.012550998556, 0.000429577973,
  125.                 -0.001801640704, 0.000353713800,
  126.                 0.000000000000, 0.000000000000,
  127.                 0.000000000000, 0.000000000000,
  128.                 0.000000000000, 0.000000000000  },
  129.  
  130.                 {    0.054415842243, 0.312871590914,
  131.                 0.675630736297, 0.585354683654,
  132.                 -0.015829105256, -0.284015542962,
  133.                 0.000472484574, 0.128747426620,
  134.                 -0.017369301002, -0.044088253931,
  135.                 0.013981027917, 0.008746094047,
  136.                 -0.004870352993, -0.000391740373,
  137.                 0.000675449406, -0.000117476784,
  138.                 0.000000000000, 0.000000000000,
  139.                 0.000000000000, 0.000000000000  },
  140.  
  141.                 {    0.038077947364, 0.243834674613,
  142.                 0.604823123690, 0.657288078051,
  143.                 0.133197385825, -0.293273783279,
  144.                 -0.096840783223, 0.148540749338,
  145.                 0.030725681479, -0.067632829061,
  146.                 0.000250947115, 0.022361662124,
  147.                 -0.004723204758, -0.004281503682,
  148.                 0.001847646883, 0.000230385764,
  149.                 -0.000251963189, 0.000039347320,
  150.                 0.000000000000, 0.000000000000  },
  151.  
  152.                 {    0.026670057901, 0.188176800078,
  153.                 0.527201188932, 0.688459039454,
  154.                 0.281172343661, -0.249846424327,
  155.                 -0.195946274377, 0.127369340336,
  156.                 0.093057364604, -0.071394147166,
  157.                 -0.029457536822, 0.033212674059,
  158.                 0.003606553567, -0.010733175483,
  159.                 0.001395351747, 0.001992405295,
  160.                 -0.000685856695, -0.000116466855,
  161.                 0.000093588670, -0.000013264203  }   };
  162.  
  163.  
  164. void main(int argc, char *argv[])
  165. {
  166.    long i;
  167.    int j;
  168.    float *data;
  169.    float *buf;
  170.    long num;
  171.    long nn;
  172.    long width;
  173.    long offset;
  174.    struct opt opt;
  175.    double ecoef[20];
  176.    double ocoef[20];
  177.    int wvshft[11][2] = {     { 0, 0 },
  178.                  { 0, 0 },
  179.                  { 0, 1 },
  180.                  { 0, 2 },
  181.                  { 0, 3 },
  182.                  { 1, 3 },
  183.                  { 1, 4 },
  184.                  { 1, 5 },
  185.                  { 1, 6 },
  186.                  { 1, 7 },
  187.                  { 1, 8 }    };
  188.  
  189.    /* load in the command line options from the application window */
  190.  
  191.    if( argc < 2 )
  192.       {
  193.       fprintf( stderr, "\nDAUBWAVE <Input file> -(options)\n" );
  194.       fprintf( stderr, "   Options:  -o <Output file> -d# -n# -t# -i# -s# -l# -h# -b# -k# -r -1 -c -u\n" );
  195.       fprintf( stderr, "   Use -? for a help file.\n" );
  196.       exit( 1 );
  197.       }
  198.    else
  199.       {
  200.       opt = InitOption( argc, argv );
  201.       }
  202.  
  203.    /* allocate the memory for storing the data and result vectors */
  204.  
  205.    data = (float *)MemAlloc( sizeof(float), MAXNUM );
  206.    buf = (float *)MemAlloc( sizeof(float), MAXNUM + 2*(opt.filter-1) );
  207.  
  208.    /* load in the data */
  209.  
  210.    num = FileLoad( data, opt );
  211.  
  212.    /* determine the largest power of two that covers the data */
  213.  
  214.    nn = SpanData( &opt.level, num );
  215.  
  216.    /* set up the coefficients for the transform */
  217.  
  218.    TranCoef( ecoef, ocoef, opt.filter );
  219.  
  220.    /* apply normalization */
  221.  
  222.    for( j=0; j<opt.filter*2; j++ )
  223.       {
  224.       ecoef[j] *= opt.normal;
  225.       ocoef[j] *= opt.normal;
  226.       }
  227.  
  228.    /* if -1 is used shift the data one to the left */
  229.  
  230.    if( (opt.shift == 1) && (opt.action != 'i') )
  231.       VectorShift( data, nn, 1, -1 );
  232.  
  233.    /*    Go through the passes of the wavelet filter and save the
  234.      high frequency component in the lower half of data and the
  235.      low frequency component in the upper half of data        */
  236.  
  237.    width = nn;
  238.    for( j=0; j<opt.level; j++ )
  239.       {
  240.  
  241.       /* avoid decomposition if only wanting to do the inverse transform */
  242.  
  243.       if( opt.action != 'i' )
  244.      {
  245.  
  246.      /* perform the wavelet transform */
  247.  
  248.      VectorCopy( buf, data, width, 2*(opt.filter-1) );
  249.      Wavelet( buf, width, 2*opt.filter, ecoef, ocoef );
  250.      VectorPermute( data, buf, width, 0 );
  251.  
  252.      /* rotate the low and high frequency numbers so they are
  253.         registered properly                     */
  254.  
  255.      if( opt.rotate )
  256.         {
  257.         VectorShift( data, width/2, wvshft[opt.filter][0], 1 );
  258.         VectorShift( data+width/2, width/2, wvshft[opt.filter][1], 1 );
  259.         }
  260.      }
  261.  
  262.       /* shift to the upper half so the decomposition can be done again */
  263.  
  264.       width /= 2;
  265.       }
  266.  
  267.    /* if only transform results are requested save and quit */
  268.  
  269.    if( opt.action == 't' )
  270.       {
  271.       FileSave( data, opt.level, nn, opt );
  272.       exit( 0 );
  273.       }
  274.  
  275.    if( opt.action == 's' )
  276.       {
  277.       FileSave( data, 0, width, opt );
  278.       offset = width;
  279.       for( j=0; j<opt.level; j++ )
  280.      {
  281.      FileSave( data + offset, j+1, offset, opt );
  282.      offset *= 2;
  283.      }
  284.       exit( 0 );
  285.       }
  286.  
  287.    /* if filtering remove either the high or low frequency component  */
  288.  
  289.    if( (opt.action == 'h') || (opt.action == 'b') )
  290.       {
  291.       for( i=0; i<width; i++ )
  292.      data[i] = 0.;
  293.       }
  294.  
  295.    if( opt.action == 'b' )
  296.       {
  297.       for( i=width*2; i<nn; i++ )
  298.      data[i] = 0.;
  299.       }
  300.  
  301.    if( opt.action == 'k' )
  302.       {
  303.       for( i=width; i<width*2; i++ )
  304.      data[i] = 0.;
  305.       }
  306.  
  307.    if( opt.action == 'l' )
  308.       {
  309.       for( i=width; i<nn; i++ )
  310.      data[i] = 0.;
  311.       }
  312.  
  313.    /* set up the inverse coefficients */
  314.  
  315.    TranCoef( ecoef, ocoef, -opt.filter );
  316.  
  317.    /* apply normalization */
  318.  
  319.    for( j=0; j<opt.filter*2; j++ )
  320.       {
  321.       ecoef[j] /= opt.normal;
  322.       ocoef[j] /= opt.normal;
  323.       }
  324.  
  325.    /* reconstruct the data set from the wavelet results            */
  326.  
  327.    for( j=0; j<opt.level; j++ )
  328.       {
  329.  
  330.       /* expand to the next largest scale to decomposite again */
  331.  
  332.       width *= 2;
  333.  
  334.       /* shift the low and high frequency components for reintegration */
  335.  
  336.       if( opt.rotate )
  337.      {
  338.      VectorShift( data, width/2, wvshft[opt.filter][0], -1 );
  339.      VectorShift( data+width/2, width/2, wvshft[opt.filter][1], -1 );
  340.      }
  341.  
  342.       /* perform the inverse wavelet transform                */
  343.  
  344.       VectorPermute( buf, data, -width, -2*(opt.filter-1) );
  345.       Wavelet( buf, width, 2*opt.filter, ecoef, ocoef );
  346.       VectorCopy( data, buf, width, 0 );
  347.       }
  348.  
  349.    /* if -1 is used shift the data one to the left */
  350.  
  351.    if( opt.shift == 1 )
  352.       VectorShift( data, nn, 1, 1 );
  353.  
  354.    /* send the results to the output */
  355.  
  356.       FileSave( data, opt.level, num, opt );
  357. }
  358.  
  359. /* VectorCopy - This routine copies the data from one vector to another.
  360.         If the padding is positive that number of values are
  361.         copied from the left side of the vector and tagged onto
  362.         the right side of the vector.  If negative then the
  363.         right values are added to the left side of the vector.    */
  364.  
  365. void VectorCopy( float *out, float *in, long length, int padding )
  366. {
  367.    long i;
  368.  
  369.    if( padding < 0 )
  370.       {
  371.       padding *= -1;
  372.  
  373.       for( i=0; i<length; i++ )
  374.      out[padding + i] = in[i];
  375.  
  376.       /* i%length is used to allow the data to be repeated multiple times
  377.      if it is shorter than padding */
  378.  
  379.       for( i=0; i<padding; i++ )
  380.      out[padding - 1 - i] = out[padding + length - 1 - i%length];
  381.       }
  382.    else
  383.       {
  384.       for( i=0; i<length; i++ )
  385.      out[i] = in[i];
  386.  
  387.       /* i%length is used to allow the data vector to be continually repeated
  388.      if it is shorter than padding */
  389.  
  390.       for( i=0; i<padding; i++ )
  391.      out[length+i] = out[i%length];
  392.       }
  393. }
  394.  
  395. /* Wavelet - This routine calculates the wavelet decomposition and
  396.          returns the values properly registered.  The high frequency
  397.          component is saved in the first half of the vector and the
  398.          low frequency is saved in the last half of the vector.    */
  399.  
  400. void Wavelet( float *buf, long length, int width, double *ecoef,
  401.                             double *ocoef )
  402. {
  403.    long i;
  404.    int j;
  405.    int times;
  406.    float even, odd;
  407.  
  408.    /* begin the loop through the data */
  409.  
  410.    for( i=0; i<length; i += 2 )
  411.       {
  412.       even = 0.;
  413.       odd = 0.;
  414.  
  415.       for( j=0; j<width; j++ )
  416.      {
  417.      even += ecoef[j]*buf[i+j];
  418.      odd += ocoef[j]*buf[i+j];
  419.      }
  420.       buf[i] = even;
  421.       buf[i+1] = odd;
  422.       }
  423. }
  424.  
  425. /* VectorPermute - This routine takes one vector and permutes it and
  426.            places it in a second vector.  Alternate values are
  427.            placed together when the length is positive and
  428.            alternate values are split up if the length is negative.
  429.            If padding is negative then the right most values are
  430.            padded to the left of the vector else the left most
  431.            values are padded to the right of the vector.    */
  432.  
  433. void VectorPermute( float *out, float *in, long length, int padding )
  434. {
  435.    long i;
  436.    long halflength;
  437.    int offset;
  438.  
  439.    if( length < 0 )
  440.       halflength = -length/2;
  441.    else
  442.       halflength = length/2;
  443.  
  444.    if( padding < 0 )
  445.       offset = -padding;
  446.    else
  447.       offset = 0;
  448.  
  449.    if( length < 0 )
  450.       {
  451.       /* place values alternately */
  452.  
  453.       length *= -1;
  454.       for( i=0; i<halflength; i++ )
  455.      {
  456.      /* copy the even index - low pass result */
  457.  
  458.      out[offset + i*2] = in[i];
  459.  
  460.      /* copy the odd index - high pass result */
  461.  
  462.      out[offset + i*2 + 1] = in[i + halflength];
  463.      }
  464.       }
  465.    else
  466.       {
  467.       /* split the alternate values up */
  468.  
  469.       for( i=0; i<halflength; i++ )
  470.      {
  471.      out[offset + i] = in[i*2];
  472.      out[offset + i + halflength] = in[i*2+1];
  473.      }
  474.       }
  475.  
  476.    if( padding < 0 )
  477.       {
  478.       /* pad the left side */
  479.  
  480.       padding *= -1;
  481.  
  482.       /* i%length is used to allow the data to be repeated multiple times
  483.      if it is shorter than padding */
  484.  
  485.       for( i=0; i<padding; i++ )
  486.      out[padding - 1 - i] = out[padding + length - 1 - i%length];
  487.       }
  488.    else
  489.       {
  490.       /* pad the right side */
  491.  
  492.       /* i%length is used to allow the data vector to be continually repeated
  493.      if it is shorter than padding */
  494.  
  495.       for( i=0; i<padding; i++ )
  496.      out[length+i] = out[i%length];
  497.       }
  498. }
  499.  
  500. /* Vec_Shift - This routine shifts a vector either right or left depending
  501.            on the value of dir.  The left-over data values are then
  502.            wrapped around to the other end.                */
  503.  
  504. void VectorShift( float *data, long length, int shift, int dir )
  505. {
  506.    float bfs[20];  /* is good up to D40 filter */
  507.    long i;
  508.  
  509.    /* if shift distance exceeds the length then reduce the shift */
  510.  
  511.    shift %= length;
  512.  
  513.    if( dir == 1 )      /* rotate upward */
  514.       {
  515.       for( i=0; i<shift; i++ )
  516.      bfs[i] = data[length-shift+i];
  517.       for( i=0; i<length-shift; i++ )
  518.      data[length-i-1] = data[length-shift-i-1];
  519.       for( i=0; i<shift; i++ )
  520.      data[i] = bfs[i];
  521.       }
  522.    if( dir == -1 )      /* rotate downward */
  523.       {
  524.       for( i=0; i<shift; i++ )
  525.      bfs[i] = data[i];
  526.       for( i=0; i<length-shift; i++ )
  527.  
  528.      data[i] = data[i+shift];
  529.       for( i=0; i<shift; i++ )
  530.      data[length-shift+i] = bfs[i];
  531.       }
  532. }
  533.  
  534. /* SpanData - This routine calculates the minimum length which is spans
  535.           the data and is also a power of 2.            */
  536.  
  537. long SpanData( int *iterates, long number )
  538. {
  539.    int loops = 0;
  540.    long nn = 1;
  541.  
  542.    while( 1 )
  543.       {
  544.       if( nn >= number )
  545.      break;
  546.       nn <<= 1;
  547.       loops++;
  548.       }
  549.  
  550.    /* the analysis can not exceed the scale of the data */
  551.  
  552.    if( *iterates == 0 )
  553.       {
  554.       *iterates = loops;
  555.       }
  556.    else
  557.       {
  558.       if( *iterates > loops )
  559.      {
  560.      fprintf( stderr, "Data not long enough for the number of iterations specified\n" );
  561.      fprintf( stderr, "Max # of iterations possible is %d.\n", loops );
  562.      exit( 1 );
  563.      }
  564.       }
  565.  
  566.    return( nn );
  567. }
  568.  
  569. /* FileOpen - This routine opens up a file and checks for errors */
  570.  
  571. FILE *FileOpen( char *file, char *opt )
  572. {
  573.    FILE *ptr;
  574.  
  575.    if( (ptr = fopen( file, opt )) == NULL )
  576.       {
  577.       fprintf( stderr, "Can not open the file %s\n", file );
  578.       exit( 1 );
  579.       }
  580.    return( ptr );
  581. }
  582.  
  583. /* File_Load - This routine loads in a data file of floating point
  584.            numbers.  If the file does not start with a floating
  585.            point number, then the file is searched for the first
  586.            line that begins with a floating point number.  The
  587.            value returned is the number of points read.        */
  588. long FileLoad( float *data, struct opt opt )
  589. {
  590.    FILE *in;
  591.    long i;
  592.    long len;
  593.    float tmp;
  594.  
  595.    /* get the input from stdin if the unix option is specified */
  596.  
  597.    if( opt.unix )
  598.       {
  599.       i = 0;
  600.       while( i < MAXNUM )
  601.      {
  602.      if( fscanf( stdin, "%g", &tmp ) == EOF )
  603.         break;
  604.      data[i++] = tmp;
  605.      }
  606.       return( i );
  607.       }
  608.  
  609.    /* open the input file */
  610.  
  611.    if( opt.compress )
  612.       {
  613.       in = FileOpen( opt.infile, "rb" );
  614.       fseek( in, 0L, 2 );
  615.       len = ftell( in );
  616.       fseek( in, 0L, 0 );
  617.       if( len%sizeof( float ) )
  618.      {
  619.      fprintf( stderr, "Error in reading input file.\n" );
  620.      exit( 1 );
  621.      }
  622.       else
  623.      {
  624.      i = len/sizeof( float );
  625.      if( i > MAXNUM )
  626.         {
  627.         fprintf( stderr, "Input file is too large.\n" );
  628.         exit( 1 );
  629.         }
  630.      else
  631.         {
  632.         if( fread( data, sizeof( float ), i, in ) != i )
  633.            {
  634.            fprintf( stderr, "Error reading input file.\n" );
  635.            exit( 1 );
  636.            }
  637.         }
  638.      }
  639.       }
  640.    else
  641.       {
  642.       in = FileOpen( opt.infile, "r" );
  643.       i = 0;
  644.       while( i < MAXNUM )
  645.      {
  646.      if( fscanf( in, "%g", &tmp ) == EOF )
  647.         break;
  648.      data[i++] = tmp;
  649.      }
  650.       }
  651.  
  652.    fclose( in );
  653.    return( i );
  654. }
  655.  
  656. /* FileSave - This routine saves the results to the output using the
  657.           appropriate extension.                    */
  658.  
  659. void FileSave( float *data, int level, long length, struct opt opt )
  660. {
  661.    long i;
  662.    FILE *out;
  663.  
  664.    /* send the output to stdout if the unix option is specified */
  665.  
  666.    if( opt.unix )
  667.       {
  668.       for( i=0; i<length; i++ )
  669.      {
  670.      if( fprintf( stdout, "%g\n", data[i] ) == EOF )
  671.         {
  672.         fprintf( stderr, "Error writing output file.\n" );
  673.         exit( 1 );
  674.         }
  675.      }
  676.       return;
  677.       }
  678.  
  679.    /* construct the output name */
  680.  
  681.    if( strlen( opt.root ) != NULL )
  682.       {
  683.       opt.extension[0] = '.';
  684.  
  685.       /* use the action type for the first byte */
  686.  
  687.       if( opt.action == 's' )
  688.      {
  689.      opt.extension[1] = HyperHex( opt.level );
  690.      }
  691.       else
  692.      {
  693.      opt.extension[1] = opt.action;
  694.      }
  695.  
  696.       /* use the wavelet type for the second byte */
  697.  
  698.       opt.extension[2] = HyperHex( 2*opt.filter );
  699.  
  700.       /* use the level of the decomposition for the third byte */
  701.  
  702.       opt.extension[3] = HyperHex( level );
  703.       opt.extension[4] = 0x00;
  704.       strcpy( opt.outfile, opt.root );
  705.       strcat( opt.outfile, opt.extension );
  706.       }
  707.  
  708.    /* open the output file and save */
  709.  
  710.    if( opt.compress )
  711.       {
  712.       out = FileOpen( opt.outfile, "wb" );
  713.       if( fwrite( data, sizeof( float ), length, out ) != length )
  714.      {
  715.      fprintf( stderr, "Error writing ouput file.\n" );
  716.      exit( 1 );
  717.      }
  718.       }
  719.    else
  720.       {
  721.       out = FileOpen( opt.outfile, "w" );
  722.       for( i=0; i<length; i++ )
  723.      {
  724.      if( fprintf( out, "%g\n", data[i] ) == EOF )
  725.         {
  726.         fprintf( stderr, "Error writing output file.\n" );
  727.         exit( 1 );
  728.         }
  729.      }
  730.       }
  731.    fclose( out );
  732. }
  733.  
  734. /* MemAlloc - This routine allocates memory, checks for errors and
  735.           sets it equal to zero.                    */
  736.  
  737. void *MemAlloc( int size, long length )
  738. {
  739.    char *ptr;
  740.    long i;
  741.  
  742.    ptr = (char *)malloc( size*length );
  743.    if( ptr == NULL )
  744.       {
  745.       fprintf( stderr, "Error in allocating memory.\n" );
  746.       exit( 1 );
  747.       }
  748.  
  749.    /* clear the memory */
  750.  
  751.    for( i=0; i<size*length; i++ )
  752.       {
  753.       ptr[i] = 0x00;
  754.       }
  755.  
  756.    return( (void *)ptr );
  757. }
  758.  
  759. /* TranCoef - This routine sets up the coefficients for the wavelet
  760.           transform.  If filt is negative then the coefficients are
  761.           set up to do the inverse transform.            */
  762.  
  763. void TranCoef( double *even, double *odd, int filter )
  764. {
  765.    int i;
  766.  
  767.    if( filter > 0 )
  768.       {
  769.       for( i=0; i<2*filter; i++ )
  770.      {
  771.      /* low pass */
  772.  
  773.      even[i] = coef[filter][i];
  774.  
  775.      /* high pass */
  776.  
  777.      odd[i] = coef[filter][2*filter - i -1];
  778.      if( i%2 )
  779.         odd[i] *= -1;
  780.      }
  781.       }
  782.    else
  783.       {
  784.       filter *= -1;
  785.       for( i=0; i<filter; i++ )
  786.      {
  787.  
  788.      /* inverse even index */
  789.  
  790.      even[i*2] = coef[filter][filter*2 - 2*(i+1)];
  791.      even[i*2+1] = coef[filter][2*i+1];
  792.  
  793.      /* inverse odd index */
  794.  
  795.      odd[i*2] = coef[filter][filter*2 - (2*i+1)];
  796.      odd[i*2+1] = -coef[filter][i*2];
  797.      }
  798.       }
  799. }
  800.  
  801. /* InitOption - This routine takes the command line options and sets the
  802.         appropriate values within the structure 'opt'        */
  803.  
  804. struct opt InitOption( int argc, char *argv[] )
  805. {
  806.    int i;
  807.    int lev;
  808.    char *ptr;
  809.    struct opt opt;
  810.  
  811.    /* set up the default values */
  812.  
  813.    strcpy( opt.infile, "noname.dat" );
  814.    strcpy( opt.outfile, "" );
  815.    strcpy( opt.root, "" );
  816.    opt.filter = 2;
  817.    opt.compress = 0;
  818.    opt.normal = 1.;
  819.    opt.action = 't';
  820.    opt.level = 0;
  821.    opt.shift = 0;
  822.    opt.unix = 0;
  823.    opt.rotate = 0;
  824.  
  825.    for( i=1; i<argc; i++ )
  826.       {
  827.       if( argv[i][0] == '-' )
  828.      {
  829.      switch( argv[i][1] )
  830.         {
  831.         case 'o':
  832.         case 'O':  strcpy( opt.outfile, argv[i+1] );
  833.                i++;
  834.                break;
  835.         case '?':  HelpFile();
  836.                exit( 0 );
  837.                break;
  838.         case '1':  opt.shift = 1;
  839.                break;
  840.         case 'n':
  841.         case 'N':  lev = atoi( &argv[i][2] );
  842.                        if( (lev < 1) || (lev > 3) )
  843.                           {
  844.                           fprintf( stderr, "%s is not a valid option.\n", argv[i] );
  845.                           exit( 1 );
  846.                           }
  847.                        if( atoi( &argv[i][2] ) != 2 )
  848.               opt.normal = 1.41421356237309504880;
  849.                if( atoi( &argv[i][2] ) > 2 )
  850.               opt.normal = 1/opt.normal;
  851.                break;
  852.         case 'u':
  853.         case 'U':  opt.unix = 1;
  854.                break;
  855.         case 'r':
  856.         case 'R':  opt.rotate = 1;
  857.                break;
  858.         case 'd':
  859.         case 'D':  opt.filter = atoi( &argv[i][2] );
  860.                if( (opt.filter%2 == 1) || (opt.filter < 2) ||
  861.                           (opt.filter > 20) )
  862.               {
  863.               fprintf( stderr, "%s is not a valid option.\n", argv[i] );
  864.               exit( 1 );
  865.               }
  866.                else
  867.               {
  868.               opt.filter /= 2;
  869.               }
  870.                break;
  871.         case 'c':
  872.         case 'C':  opt.compress = 1;
  873.                break;
  874.         case 'l':
  875.         case 'L':  opt.action = 'l';
  876.                lev = atoi( &argv[i][2] );
  877.                if( lev > 0 )
  878.               opt.level = lev;
  879.                break;
  880.         case 'h':
  881.         case 'H':  opt.action = 'h';
  882.                lev = atoi( &argv[i][2] );
  883.                if( lev > 0 )
  884.               opt.level = lev;
  885.                break;
  886.         case 'b':
  887.         case 'B':  opt.action = 'b';
  888.                lev = atoi( &argv[i][2] );
  889.                if( lev > 0 )
  890.               opt.level = lev;
  891.                break;
  892.         case 'k':
  893.         case 'K':  opt.action = 'k';
  894.                lev = atoi( &argv[i][2] );
  895.                if( lev > 0 )
  896.               opt.level = lev;
  897.                break;
  898.         case 'i':
  899.         case 'I':  opt.action = 'i';
  900.                lev = atoi( &argv[i][2] );
  901.                if( lev > 0 )
  902.               opt.level = lev;
  903.                break;
  904.         case 't':
  905.         case 'T':  opt.action = 't';
  906.                lev = atoi( &argv[i][2] );
  907.                if( lev > 0 )
  908.               opt.level = lev;
  909.                break;
  910.         case 's':
  911.         case 'S':  opt.action = 's';
  912.                lev = atoi( &argv[i][2] );
  913.                if( lev > 0 )
  914.               opt.level = lev;
  915.                break;
  916.         default : fprintf( stderr, "%s not a valid option.\n", argv[i] );
  917.               exit( 1 );
  918.         }
  919.      }
  920.       else
  921.      {
  922.      strcpy( opt.infile, argv[i] );
  923.      }
  924.       }
  925.  
  926.    /* set up the root and extension of the output files */
  927.  
  928.    if( strlen( opt.outfile ) == 0 )
  929.       {
  930.       strcpy( opt.outfile, opt.infile );
  931.       ptr = strchr( opt.outfile, '.' );
  932.       if( ptr )
  933.      ptr[0] = 0x00;
  934.       strcpy( opt.root, opt.outfile );
  935.       }
  936.    else
  937.       {
  938.       if( strchr( opt.outfile, '.' ) == NULL )
  939.      {
  940.      strcpy( opt.root, opt.outfile );
  941.      }
  942.       else
  943.      {
  944.      if( opt.action == 's' )
  945.         {
  946.         fprintf( stderr, "Can not specify the output file extension with \
  947.              option -s\n" );
  948.         }
  949.      }
  950.       }
  951.  
  952.    return( opt );
  953. }
  954.  
  955. /* HyperHex - This routine takes an integer between 0 and 36 and converts
  956.           it into a character that goes from 0-9,A-Z.        */
  957.  
  958. char HyperHex( int value )
  959. {
  960.    if( (value > 35) || (value < 0) )
  961.       {
  962.       fprintf( stderr, "Value too big to specify the file extension properly.\n" );
  963.       exit( 1 );
  964.       }
  965.    if( value < 10 )
  966.       {
  967.       return( '0' + value );
  968.       }
  969.    else
  970.       {
  971.       return( 'A' + value - 10 );
  972.       }
  973. }
  974.  
  975. /* HelpFile - This routine prints out a help file to the screen     */
  976.  
  977. void HelpFile( void )
  978. {
  979.    printf( "DAUBWAVE - General purpose orthogonal wavelet program.\n" );
  980.    printf( "Author: Steven Gollmer (nls@mace.cc.purdue.edu)   Date: October 23, 1992\n\n" );
  981.    printf( "DAUBWAVE <input> (options)\t Default Options: -d4 -n2 -t\n" );
  982.    printf( "Options:\n" );
  983.    printf( "-o <output>\n\tSpecify the root or the whole name of the output file.\n" );
  984.    printf( "-d#\tSpecify the wavelet transform d2 through d20.\n" );
  985.    printf( "-n#\t1/2 Normalization. 1-On inverse, 2-Root on both, 3-On transform.\n" );
  986.    printf( "-t#\tPerform the transform to # levels of the hierarchy.\n" );
  987.    printf( "-i#\tPerform the inverse transform to # levels of the hierarchy.\n" );
  988.    printf( "-s#\tSame as -t# but save each hierarchy as a separate file.\n" );
  989.    printf( "-l#\tPerform a low pass filter on the data transformed to the # level.\n" );
  990.    printf( "-h#\tPerform a high pass filter on the data transformed to the # level.\n" );
  991.    printf( "-b#\tPerform a band pass filter on the data transformed to the # level.\n" );
  992.    printf( "-k#\tPerform a notch filter on the data transformed to the # level.\n" );
  993.    printf( "\t(Analyze to the highest level possible if # is not specified.)\n" );
  994.    printf( "-r\tCorrelate coefficients at different levels by shifting results.\n" );
  995.    printf( "-1\tShift the data by one before the transform is performed.\n" );
  996.    printf( "-c\tUse IEEE binary floating point for input and output.\n" );
  997.    printf( "-u\tImplement Unix type redirection for input and output.\n" );
  998.    printf( "Reference:\tPress, William H., 'Numerical Recipes for Fortran, 2nd Ed.'\n" );
  999.    printf( "\t\tNew York: Cambridge University Press, 1992.\n" );
  1000. }
  1001.